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,-rt I This paper gives a summary of basic concepts of density-functional theory (DFT) 

^1 ' and its use in state-of-the-art computations of complex processes in condensed matter 

physics and materials science. In particular we discuss how microscopic growth param- 
eters can be determined by DFT and how on this basis macroscopic phenomena can 
be described. To reach the time and length scales of realistic growth conditions, DFT 

^ I results are complemented by kinetic Monte Carlo simulations. 
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0^ . INTRODUCTION 

The microscopic processes governing epitaxial growth typically are rather complex 

_ since they involve the making and breaking of chemical bonds as well as the dynamics 

rH ■ of the nuclei. With recent progress in the developments of new methods and techniques 

O ■ and the availability of faster computers, density-functional theory (DFT) calculations 

have evolved into a powerful tool to study growth phenomena (as well as other complex 

processes in condensed matter physics, materials science, and chemistry). The rate of 

r> I a microscopic process j that occurs during growth, such as diffusion, usually has the 

c^ ■ form r^-'^ = Tq exp{—E^ /k-BT), where Fg is the effective attempt frequency, T the 

temperature, /cb the Boltzmann constant, and E^ is the energy barrier that needs 

to be overcome for the event j to take place. This equation reflects the idea that 

an adatom experiences many stable and metastable sites at the surface, and that the 

diffusive motion brings it from one minimum to an adjacent one on the free energy 

surface in the space of the reaction coordinates. The effective attempt frequency Fq 

contains the terms due to the adatom and substrate vibrations (see. Ref. 1 for more 

details). 

It is at the heart of theoretical studies of surface diffusion and growth phenomena 
to calculate the ground-state total energy of the adsorbate system for a dense mesh of 
adatom positions. This yields the so-called potential-energy surface (PES) which is the 



potential energy experienced by the diffusing adatom, 

E''^''{X^,Y^)= niin E*°*(Xad,yad,^ad,{R/}) , (1) 

where i?*°*(Xad, ^ad, ^ad, {R-/}) is the ground-state energy of the many-electron system 
(also referred as the total energy) at the atomic configuration (Xad, ^ad, -^ad, {R-/})- 
According to Eq. (|l|) the PES is the minimum of the total energy with respect to the 
2;-coordinate of the adatom and all coordinates of the substrate atoms {R/}. Assum- 
ing that vibrational effects can be neglected, the minima of the PES represent stable 
and metastable sites of the adatom. Note, this PES refers to slow motion of nuclei 
and assumes that for any atomic configuration the electrons are in their respective 
ground state. Thus, it is assumed that the dynamics of the electrons and of the nuclei 
are decoupled. This is the Born-Oppenheimer approximation which for not too high 
temperatures is usually well justified. 

Now consider all possible paths j to get from one stable or metastable adsorption 
site, Rad, to an adjacent one, Rad'- The energy difference E^^ between the energy 
at the saddle point along j and the energy of the start geometry is the barrier for 
this particular path. The diffusion barrier then is the minimum value of all E^^ of all 
possible paths which connect Rad and Rad', and the lowest energy saddle point is called 
the transition state. We note this definition strictly applies only to cases where the 
vibrational energy is negligible, which is typically justified when the diffusion barrier 
is not too small, and the temperature is not too high. We also note that for non zero 
temperatures, i.e., when vibrations are thermally excited, vibrational entropy needs to 
be considered, which will enter the attempt frequency of the hopping rate. Although 
often only the path with the most favorable energy barrier is important, it may happen 
that several paths exist with comparable barriers or that the PES consists of more than 
one sheet (e.g., Ref. 2). Then the effective barrier measured in an experiment (or a 
molecular dynamics simulation) represents a proper average over all possible pathways. 
The above description obviously applies to simple jumps of an adatom (i.e., diffusion 
by hopping), but it also holds for the diffusion by atomic exchange where the diffusing 
adatom displaces a substrate atom.^"^ 

Aiming at a calculation and understanding of the PES of a diffusing atom, it is 
obvious that the interplay between the breaking and making of chemical bonds and 
the atomic relaxations needs to be accounted for by an accurate, quantum-mechanical 
description of the many-electron system. This can be achieved by modern density- 
functional theory calculations that combine electronic self-consistency and efficient ge- 
ometry optimization (see for example Ref. 8-10 and references therein). 

In the following two sections we describe briefly the basic concepts of density- 
functional theory (DFT) and of its application to surface problems. We then discuss 
some examples where DFT has been used to calculate growth parameters and compare 
the results with experiments. Finally, an example is given that shows how DFT cal- 
culations can be used in combination with kinetic Monte Carlo simulations to describe 
and analyze epitaxial growth. 

DENSITY-FUNCTIONAL THEORY: BASIC CONCEPTS 

The total energy of an A^-electron, poly-atomic system is given by the expecta- 
tion value of the many-particle Hamiltonian using the many-body wave-function of the 



electronic ground state. For a solid or a surface the calculation of such an expecta- 
tion value is impossible when using a wave-function approach. However, as has been 
shown by Hohenberg and Kohn/^ the ground-state total energy can also be obtained 
without explicit knowledge of the many-electron wave-function, but by minimizing an 
energy functional E[n\. This is the essence of density-functional theory (DFT), which is 
primarily (though in principle not exclusively) a theory of the electronic ground state, 
couched in terms of the electron density n(r) instead of the many-electron wave-function 

The important theorem of Hohenberg and Kohn^^ (see also Levy^^) tells: The spec- 
ification of a ground state density n(r) determines the corresponding external potential 
f^^*(r) uniquely (to within an additive constant), 

n(r) ^t;""*(r) . (2) 

The external potential v^^^{y) is typically (and definitely for our purpose here) the 
Coulomb potential due to the nuclei. While the other direction [v^^^{v) -^ n(r)] is well 
known to exist, because f'^^*(r) determines the many-particle Hamiltonian, Eq. (^ is 
less obvious. In other words, Hohenberg and Kohn realized that for the ground state 
the functional n{v) = n[$] = (^E'|I]j5(r — ri)|\E') can be inverted, i.e., \E' = \l'[n(r)]. 
With the help of this theorem the variational problem of the many-particle Schrodinger 
equation transforms to a variational problem of an energy functional: 

Eo < i^^lHl^) = E4^[n]] = E4n] . (3) 

Here Eq is the energy of the ground state, and Ejj[n\ = J drv^^^{r)n{r) + G[n]. In 
this functional n{r) is the variable (the electron ground-state density of any A^-electron 
system), and f''^*(r) is kept fixed. G[n] is a universal functional independent of the 
system, i.e., independent oiv'^^^{r). For example, G[n\ is the same for an H-atom, a CO- 
molecule, a solid etc. The main advantage of this approach is that n{r) only depends 
on three variables while \E'({rj}) depends on many variables (the 3A^ coordinates of all 
electrons) . ^^ Thus, it is plausible that the variational problem of Ey [n] is easier to solve 
than that of {'^\H\'^), yet the result for the ground-state energy and the ground state 
electron density will be the same. The total energy entering Eq. (p is^ 



.14 



E'°\{Rj}) = Eo{{Rj}) + 1 E iB^'^B I ' (4) 

where {Rj} includes all atoms, and Zj is the nuclear charge. 

An important problem remains, namely that an explicit form of the functional 
G[n] is unknown. Earlier work (in particular the Thomas- Fermi approach) had shown 
that the treatment of the kinetic energy (^| — |V^|^) is of particular importance and 
Kohn and Sham^^ therefore wrote the energy functional in the form 

Ey[n] = Ts[n] + /"rfr t;""*(r)r2(r) + - f drv^{r)n{r) + E^'^ln] , (5) 

where Ts[n] is the functional of the kinetic energy of a system of non- interacting 
electrons with density n(v), and f^(r) = / dr' ll_J, is the Hartree potential that de- 
scribes the electrostatic interaction between electrons. ii'^'^[n] is the so-called exchange- 
correlation functional. It accounts for the Pauli principle, dynamical correlations due 



to the Coulomb repulsion, and the correction of the self-interaction included for con- 
venience in the Hartree term. With Eq. (^) the problem of the unknown functional 
G[n\ is transformed to one involving Ts[n\ and E^^[n\. We note in passing that the 
functional defined by Eq. (|]) can be also modified by adding terms which vanish at 
the correct electron density. Such a functional E^[n] may converge faster towards the 
ground state or may depend less sensitive on the input density. The latter implies that 
the input density does not need to be very good, yet the resulting energy represents 
an acceptable approximation for the correct total energy (see e.g., Ref. 16). Although 
the functional Ts[n\ is not known explicitly in a mathematically closed form, it can be 
evaluated exactly by using the following "detour" proposed by Kohn and Sham. The 
variational principle applied to Eq. (H) leads to 

where /i is the Lagrange multiplier associated with the requirement of a constant particle 
number and thus equals the electron chemical potential. The effective potential is 
defined as 
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with f^^(r) = 5E^'^[n\/5n{Y), and n{v) is a ground-state density of any non-interacting 
electron system, i.e., 

n{v)=Y.U\U^)? , (9) 

where we introduced the occupation numbers /j. Because Ts[n] is the kinetic energy 
functional of non- interacting electrons, Eq. (J^) (together with Eq. (^) is solved by 



.lv= + ,. 



e,0,(r) . (10) 



These are the Kohn-Sham equations, that are to be solved self-consistently together 
with Eqs. (|) and (||). In principle this gives the exact ground-state electron density and 
total energy of a system of interacting electrons. However, the functional i?^'^[n] is still 
unknown. Some general properties of this functional and values for some special cases 
are known. Detailed and very accurate understanding exists for systems of constant 
electron density. The asymptotic behavior at low and high densities is given by expres- 
sions derived by Wigner^^ and Gell-Mann and Brueckner^^ and for intermediate densi- 
ties quantum Monte Carlo calculations have been performed by Ceperley and Alder. ^^ 
This gives the simple curve shown in Fig. |l], and this result for e'"^(rs) := ^ij^xi''^) is 
then used in the functional 



^LDAM = /rfrn(r)6-A(^(r)) , (11^ 



which is the local-density approximation (LDA).^^ Thus, in the LDA the many-body 
effects are included such that for a homogeneous electron gas the treatment is exact 
and for an inhomogeneous system exchange and correlation are treated by assuming 
that the system can be composed from many small systems with a locally constant 
density. 



-10.0 



>. 

DC 

u 

X 



-20.0 



-30.0 




100.0 



rs(bohr) 



Figure 1. Exchange-correlation energy per particle, e^'^, of homogeneous electron gases with density 
parameters Tg. The electron density and the density parameter are related by n = t'^'''^- 

The LDA can be improved by including the dependence on the density gradient 
which leads to the generalized gradient approximation (GGA). Several different GGA's 
were proposed in the literature ^''"^^ and have been used successfully for DFT calcu- 
lations of atoms, molecules, bulk solids, and surfaces (an overview can be found in 
Refs. 25 and 26), but also limitations have been identified for example by Mitas et 
al?^ and Umrigar and coworkers. ^^ It is by now clear that the lattice constants cal- 
culated with a GGA are typically larger than those obtained with the LDA, with the 
experimental values usually being in between. Binding energies (or cohesive energies) 
of molecules and solids are clearly improved by the GGA and energy barriers of chem- 
ical reactions are improved as well (see Ref. 29 and 30 and references therein). Still, 
for surface diffusion DFT-LDA calculations give energy barriers in good agreement 
with those deduced from experiments and with GGA calculations. Although the to- 
tal energies are changed when going from the LDA to the GGA, the changes in energy 
barriers, i.e., in total energy differences, are typically less pronounced (see e.g., Ref. 31). 



IMPLEMENTATION OF DENSITY-FUNCTIONAL THEORY TO COM- 
PUTE MICROSCOPIC PARAMETERS 



Typically there are only a few candidates for adsorption sites and possible channels 
for diffusion. This is illustrated in Fig. |^ for the fee (111) and fee (100) surfaces. For 
adatoms which are chemically similar (or equivalent) to those of the substrate the 
stable sites are those with high coordination and for hopping diffusion the transition 
state is at the bridge site. The relevant information about the PES then is obtained by 
calculating the total energy of the system with the adatom placed in those positions. 
In general, more care is necessary because the bridge site can also be a local minimum 
of the PES and the energy barrier could be in between the high coordination and the 
bridge sites. Furthermore, it is possible that the diffusion does not proceed by hopping 
but by atomic exchange. '^"'^ 



In the bulk crystal the three-dimensional periodicity can be exploited by using 
Bloch's theorem. Unfortunately, the presence of a surface and an adatom on top of it 
breaks all periodicity. The (in principle) best approach to treat such difficult situation 
is given by the Green-function method. ^^'^^ A popular approximation for (at least in 
the past) for adsorbate systems is the cluster approach. ^^ The presently most efficient 
and practical approach that was also used in the results discussed below is the supercell 
approach. The supercell may be also called a big cluster, but in contrast to conventional 
cluster calculations the supercell is periodically repeated. As a consequence, the cluster 
boundary is treated physically very accurately, and by utilizing the periodicity, i.e., the 
Bloch theorem, it is possible to use very big cells. The idea of an adatom on top of a 
substrate in the supercell approach is sketched in Fig. |^. The adatom is placed on top 
of a slab of a certain number of layers. The number of layers iVgiab must be sufficiently 
large so the adatom does not "feel" the presence of another surface on the other side of 
the slab (or at least that the quantity to be computed, e.g., a diffusion barrier, is not 
affected by the other surface). Alternatively, one could also place an adatom on either 
side of the slab; in this case, there are more symmetries in the geometry but more 
layers are needed in the slab to screen the mutual interaction between the two adatoms 
through the slab. The adequate number of layers A^siab depends on the properties that 
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Figure 2. Top view at a fee (111) and (100) surfaec. Tlic adsorption sites labeled fee, hep, and 
hollow site usually eorrespond to tlie most stable binding sites while the bridge-site is the transition 
state of a hopping diffusion. 
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Figure 3. Sketch of a supercell describing an "isolated" adatom at a surface (side view). 



one wants to calculate and the surface orientation, and careful tests must be carried 
out. For example, for the Ag(lll) surface four layers are sufficient when the adatom 
is placed on only one side of the slab, while for Al(lOO) seven layers are necessary. 

As illustrated in Fig. § the geometry repeats periodically in vertical and lateral 
directions. The lateral periodicity implies that a single adatom placed on a substrate is 
not at all a single adatom; if the cell size parallel to the surface is chosen, for example, as 
(2 X 2) we actually calculate a system with a coverage of 25 %. It is therefore important 
to test that the interaction with the neighboring adatoms can be neglected. On a (111) 
surface a cell size of (2 x 2) is usually sufficient, but sometimes larger cells [(3 x 3) 
or even (4 x 4)] are necessary. To model a diffusion event along or across a step one 
either chooses a small island on top of a substrate or a vicinal surface. The latter has 
the advantage that only one step exists in the unit cell so a smaller cell size is required 
to attain a negligible step-step interaction. The system also repeats in the vertical 
direction separated by a vacuum region. The thickness of the vacuum region must be 
tested as well, but the computational cost of a thicker vacuum region is relatively small 
compared to a larger cell size or a higher number of slab layers (for a deeper discussion 
of the above technicalities, see e.g., Ref. 35). 

Core electrons typically do not take part directly in the binding process of atoms 
in molecules and solids, and the nature of the chemical bond is mainly determined by 



the valence electrons. This is exploited by the frozen core approximation where the core 
electrons are effectively combined with the nuclei to form frozen (i.e. unpolarizable) 
ions. Still, not just the electrostatic potential but also the quantum nature of the core 
electrons is felt by the valence electrons. For example, different wave functions must be 
orthogonal and therefore the valence wave functions have nodes and oscillate in the core 
region. For practical calculations one needs to expand the wave function in a suitable 
basis and we choose a plane wave basis set^^ 

|</.,(k,r))=5:c,-k(G)|G + k) . (12) 

G 

A plane wave description of wave functions that have nodes and oscillate requires 
a very large number of plane waves. This inconvenience is cured efficiently by the 
pseudopotential approach. Modern ab initio pseudopotentials reproduce the potential 
of an atom exactly outside the core region defined by a radius r^ and are rather smooth 
inside the core region. An important requirement on a "good" pseudopotential is 
that it is transferable, which means that the pseudopotential should behave like the 
all-electron potential in a variety of different chemical situations. Pseudopotentials 
that reproduce the same charge inside the core region as the all-electron potential, 
and therefore have the same scattering properties, are referred to as norm-conserving. 
Those that are often used have been developed by Bachelet, Hamann, and Schliiter^^, 
TrouUier and Martins^^, and Gonze, Stumpf, and SchefHer.^^ Recently, Vanderbilt ^° 
proposed ab initio pseudopotentials that drop the condition of norm-conservation and 
therefore can be used with a lower number of plane waves. The gain in computer 
time due to the smaller basis set is partially compensated by the costs to calculate the 
correction required by the neglected norm-conservation. 

The electron density is calculated according to Eq. (^ as 

^(r) = EE^k/(6,(k))|0,(k,r)r (13) 

k j 

where the integration over the Brillouin zone has been replaced by a sum over a mesh 
of k-points with u^ the k-point's weight. Such replacement typically works quite ef- 
ficiently because the electron wave functions vary rather smoothly with k, and the 
main dependence usually is in the phase-factor. Thus, a certain region in k-space is 
well represented by only one k-point. A convenient scheme to construct an appropri- 
ate k-point mesh is described by Monkhorst and Pack.^^ In ab initio pseudopotential 
calculations some matrix elements and some integrals are efficiently evaluated in real 
space whereas others are efficiently evaluated in reciprocal space. The technique of fast 
Fourier transformation enables a numerically fast change from one representation to 
the other. Technical details of the computational procedures are described for example 
in Refs. 8 and 10. 



RESULTS FOR MICROSCOPIC PARAMETERS 

We now discuss some results obtained by DFT calculations that provide a deeper 
insight into the microscopic mechanisms behind growth phenomena. The main objec- 
tive is to identify the nature and to determine the energetics of diffusion processes. For 
the homoepitaxial growth of metallic systems such as Al/Al (111),^^'^^ Al/Al (100),^^''*^ 
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and Ag/Ag (100)^'^^ comprehensive DFT studies have determined diffusion barriers for 
diffusion on the ffat surface and along and across step edges with the hopping and 
the exchange mechanism. In this section the properties of Ag/Ag (111) are discussed 
and the next section extends earher DFT calculations of Al/Al (111)^^''^^ and analyzes 
actual growth conditions of mesoscopic islands on long time scales. 

Growth of one material on a different material is of particular interest for a number 
of technological applications. In such a heteroepitaxial system with usually different 
lattice constants the material to be deposited is under the influence of epitaxial strain. 
Growth of Ag on Pt (111) and Ag on a thin Ag film on Pt (111) has been the focus of a 
number of recent studies, ^^"^^ and with a lattice mismatch of 4.2 % it provides important 
information about the effects of strain during growth. Here, we will particularly discuss 
how strain affects the surface diffusion barrier. 

Only few theoretical studies of the influence of lattice mismatch on the diffusion 
barrier are present in the literature. For a metallic system we are only aware of results 
for Ag on Ag (111) where the authors of Ref. 44 find in a semi-empirical effective medium 
theory (EMT) calculation that the diffusion barrier increases under tensile strain and 
decreases under compressive strain. 

Here, we present first-principles calculations (more details are given in Ref. 47) 
studying the dependence of the diffusion barrier on the lattice constant for Ag on 
Ag(lll).^^ In the range of ±5% strain the DFT results exhibit an approximately 
linear dependence with a slope of ~ 0.7 eV (see Fig. H). The calculated diffusion 
barrier for the unstrained system, E^^~ ^ = 81 meV, is in good agreement (within 
the error margins of the experiment and the calculations) with the scanning tunneling 
microscopy (STM) results of E^^~ ^ = 97 meV. The accordance between experiment 
and theory extends to the system Ag/Pt (111) and Ag/IML Ag/Pt (111). These results 
are summarized in Table ^ In Fig. ^the DFT-LDA results are compared to those of an 
EMT study. ^^ The EMT results exhibit a linear dependence only for very small values 
of strain (±2 %) and the diffusion barrier starts to decrease for values of misfit larger 
than 3 %. Indeed, it is plausible that a decrease of the diffusion barrier occurs when 
the atoms are separated far enough that eventually bonds are broken. However, as our 
DFT-LDA results show, for Ag/Ag (111) this happens at values for the misfit that are 



120 



> 
E 



m 
ffl 

c 
_o 
"« 

3 



100 




0.95 1.00 

Relative Lattice Constant 



1.05 



Figure 4. Diffusion barrier (in meV) for Ag on Ag (111) as function of strain. The circles are 
DFT-LDA results from Ref. 47 and the squares are EMT results from Ref. 44. 
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larger than 5 %. Additionally, when comparison with experiment is possible [i.e., Ag on 
Ag (111), and Ag on a monolayer Ag on Pt (111)] the EMT results are off by a factor 
that varies from 1.2 to 2. 

The DFT results in Fig. ^ were obtained with the LDA for the exchange-correlation 
functional and test calculations show that GGA increases the diffusion barrier by no 
more than 5 — 10 %. This has also been found for Pt on Pt (111) and Ag on Ag (100) 
(Refs. 48 and 31, respectively). The general trend of an increasing energy barrier 
for hopping diffusion with increasing lattice constant is quite plausible (for exchange 
diffusion, see Ref. 7). Smaller lattice constants correspond to a reduced corrugation of 
the surface, and as a result of a large compression the atom is not bonded much stronger 
at the adsorption sites than at the bridge site. In contrast, when the surface is stretched 
the corrugation increases and the adsorption energy at the three-fold coordinated hollow 
sites increases. This picture will change when the strain is so large that bonds are broken 
and then it is expected that the hopping diffusion barrier will start to decrease again 
at very large tensile strain. 

It is worth noting that the diffusion barrier for Ag on top of a pseudomorphic 
layer of Ag on Pt (111) is substantially lower than that for Ag on Ag (111). A question 
that arises is whether this reduced diffusion barrier is a result of the compressive strain 
or electronic effects due to the Pt substrate. The diffusion barrier for Ag on Ag(lll) 
with a lattice constant that is compressed to the value of the lattice constant for Pt 
is E^^~ ^ = 60 meV while that for Ag on Pt (111) (also with the Pt lattice constant 
of 3.92 A obtained from DFT) is E^^~ = 65 meV. The agreement of these two 

values suggests that the reduction of the diffusion barrier for Ag on a layer of Ag on 
Pt (111) is mainly a strain effect and that the diffusion barrier on top of a layer of Ag 
is essentially independent of the substrate underneath. 

We note in passing that the increase of the hopping diffusion barrier with tensile 
strain is also to be expected (and found) for the more open (100) surface. On the other 
hand, for exchange diffusion the slope of the energy barrier as a function of strain has 
the opposite sign and it has been argued^ that this behavior and the large surface stress 
at late 5d transition metals actuate exchange diffusion experimentally found for Ir (100) 
and Pt(lOO). 

AB INITIO KINETIC MONTE CARLO SIMULATIONS 

The time between two successful diffusion events is often of the order of nanosec- 
onds. Since molecular dynamics (MD) calculates all unsuccessful attempts (usually 
~ 10^) explicitly, a typical MD simulation can cover at most times of some picosec- 
onds, possibly some nanoseconds. Therefore, although MD simulations can provide 
important insight into elementary microscopic mechanisms, they normaly cannot be 



System 


Experiment (Ref. 44) 


EMT (Ref. 44) 


DFT (Ref. 47) 


Ag/Pt(lll) 
Ag/IML Ag/Pt(lll) 
Ag/Ag(lll) 


157 

60 

97 


81 
50 
67 


150 

65 

81 



Table 1. Diffusion barriers (in meV) for Ag on Pt (111), Ag on one monolayer (ML) 
Ag on Pt (111), and Ag on Ag (111). 
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used for growth studies. Instead, the method of choice for studying the spatial and 
temporal development of growth is kinetic Monte Carlo (KMC). The key idea behind 
KMC is to describe stochastic processes (such as deposition, diffusion, desorption, etc.) 
on the microscopic scales by rates and thus to avoid the explicit calculation of unsuc- 
cessful attempts. Yet, the result of a KMC study will be the same as that of an MD 
simulation, provided that the underlying PES is the same. The strategy of KMC can 
be summarized as follows: 

1) Determination of all processes j that possibly could take place with the actual 
configuration of the system. 

2) Calculation of the total rate R = Y.j r*-"'\ where the sum runs over the possible 
processes [see step 1)]. Deposition is accounted for in this description by the 
deposition rate F . 

3) Choose a random number p in the range (0, 1]. 

4) Find the integer number / for which 

i-i I 

^r(^) <pi?<^r(^) . (14) 

5) Let process / to be taken place. 

6) Update the simulation time t := t + At with At = —ln{p)/R . 

7) Go back to step 1). 

KMC simulations have been used to study crystal growth of semiconductors (e.g., 
Refs. 49-51) and metals (e.g., Refs. 52-55). However, most of these studies have been 
based on restrictive approximations. For example, the input parameters, such as acti- 
vation barriers, have been treated as effective parameters determined rather indirectly, 
e.g. from the fitting of experimental quantities, like intensity oscillations in helium 
atom scattering (HAS) measurements, in reflection high energy electron diffraction 
(RHEED) experiments, or island densities from STM pictures. Thus, the connection 
between these parameters and the microscopic nature of the processes may be somewhat 
uncertain. Often even the surface structure was treated incorrectly, i.e., the simulation 
was done on a simple cubic lattice while the system of interest had an fee or bcc struc- 
ture. Despite these approximations, such studies have provided significant qualitative 
and in some cases also quantitative insight into growth phenomena. The next better 
approach is to use semi-empirical calculations such as the embedded atom method or 
effective medium theory to evaluate the PES for KMC simulations of growth. ^^"^^ The 
best, but also most elaborate approach to obtain the PES was described in the previ- 
ous sections. In the following the DFT results for Al on Al (111) obtained by Stumpf 
and Scheffler^^''^^ are utilized for KMC simulations. Thus, it is our aim to perform a 
realistic simulation that takes into account the correct structure of the system and rate 
constants determined from accurate ab initio calculations. 

On the (111) surface of an fee crystal there are two different types of close-packed 
steps, shown in Fig. ||. They are labeled according to their {100} and {111} facets, 
referring to the plane passing through the adatom of the step and the neighbor atom 
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Figure 5. Side view (upper panel) and top view (lower panel) of the two types of close-packed 
steps on the (111) surface of an fee crystal. 



of the substrate (often these two steps are labeled A and B, respectively). Experi- 
mentally it has been shown that for Pt (111)^^ and Ir(lll)^*' these two steps behave 
quite differently with respect to surface diffusion and growth. For Al(lll) the DFT 
calculations^^'^^ predict that the formation energies of the two steps are different with a 
lower energy cost for the formation of the {111} faceted step than of the {100} faceted 
step: 0.232 eV per atom vs. 0.248 eV per atom. This difference affects the equilibrium 
shape of the islands as determined by the Wulff construction. Because more open steps 
have a higher formation energy, one expects in thermodynamic equilibrium and at not 
too high temperatures hexagonally shaped islands where the edges alternate between 
those with a {100} and a {111} microfacet, the latter being longer. 

We now like to analyze typical growth conditions, i.e., a situation far from equilib- 
rium where kinetic processes are dominant. For Al on Al (111) Stumpf and SchefHer'^^'^^ 
analyzed microscopic diffusion processes and in particular determined the activation en- 
ergies Ed for: 

(i) diffusion of a single adatom on the flat surface: E^ = 0.04 eV; 

(ii) diffusion from upper to lower terraces which was found to proceed by an exchange 
with a step-edge atom: E^ = 0.06 eV for the {100} step and E^ = 0.08 eV for 
the {111} step; 

{Hi) diffusion parallel to the {100} step via hopping: E^ = 0.32 eV (0.44 eV for 
exchange) ; 

diffusion parallel to the {111} step via exchange: E^ = 0.42 eV (0.48 eV for 
hopping) . 
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Figure 6. A surface area of (1718 x 1488) A^ (half of the simulation array) at four different 
substrate temperatures. The deposition rate was 0.08 ML/s and the coverage in each picture is 6 = 
0.04 ML. 

The DFT calculations^^ give that the binding energy of two adatoms in a dimer is 
0.58 eV, and we therefore assume that dimers, once they are formed, are stable, i.e., they 
will not dissociate. Moreover, in the lack of reliable information we assume that dimers 
are immobile. We note that the reported value^^'^^ for the self-diffusion energy barrier 
is rather low (0.04 eV) and comparable to the energy of optical phonons of Al(lll) 
(0.03 — 0.04 eV).^^ Thus, simulations at room temperature may not be reliable because 
the concept of single jumps between nearest neighbor sites is no more valid. A single 
optical phonon can furnish enough energy to an adatom for leaving its adsorption site 
and diffusing on the fiat surface. At room temperature the level population of optical 
phonon is high and the adatoms have practically no saddle point and migrate freely on 
the fiat surface. We therefore limited our study to substrate temperatures T < 250 K. 

We adopt periodic boundary conditions, and our rectangular simulation area is 
compatible with the geometry of an fee (111) surface. The dimensions of the simulation 
area are 1718 x 2976 A^. These dimensions are a critical parameter and it is impor- 
tant to ensure that the simulation area is large enough that artificial correlations of 
neighboring cells do not affect the formation of growth patterns. The mean free path 
A of a diffusing adatom before it meets another adatom with possible formation of a 
nucleation center or is captured by existing islands should be smaller than the linear 
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dimension of the simulation array. Since A is proportional to (D/FY^^ (Ref. 62), we 
have that (with F = 0.08 ML/s) A ~ 50 A for T = 50 and gets as large as ~ 10^ 
A for T = 250 K. We see that our cell is large enough (for the imposed deposition 
rate) for T < 150K, whereas at higher temperatures the dimensions of the cell are 
too small, i.e., for T > 150K the island density is determined by the simulation array 
rather than the physics. Nevertheless, the island shape is determined by local processes 
(edge diffusions) and is still meaningful. 

In the KMC program two additional insights extracted from the DFT calculations 
are included: (i) the attractive interaction between steps and single adatoms, and (u) 
the fact that diffusion processes take place via different mechanisms (hopping or ex- 
change). Particularly the second point plays an important role in our investigation. 
In several KMC simulations of epitaxial growth the attempt frequency of the diffusion 
rate has the same value for all the processes, and this value lies usually in the range of 
a typical optical phonon vibration or the Debye frequency. However, this assumption 
may not be right. First, processes with larger activation barriers may have a larger 
attempt frequency than processes with smaller energy barriers. This is a consequence 
of the compensation effect described, e.g., in Ref. 63. Moreover, processes as hopping 
and exchange that involve a different number of particles and different bonding con- 
figurations may also be characterized by different attempt frequencies. This has been 
observed ^^"^^ for several systems (Rh, Ir, Pt) and implies that the attempt frequency 
for exchange diffusion can be larger by up to two orders of magnitude than that for 
hopping. For Al surfaces calculations with the embedded atom method^^ showed a 
difference of pref actors of one order of magnitude. 

The results of the ab initio KMC simulations are shown in Fig. ^ for coverages 
of 6 = 0.04 ML. When the substrate temperature is 50 K during growth the shape 
of the islands is highly irregular and indeed fractal. Adatoms which reach a step 
cannot leave it anymore and they even cannot diffuse along the steps. Thus, at this 
temperature ramification takes place into random directions, and island formation can 
be understood in terms of the so-called hit and stick model (see also Ref. 69). At a 
growth temperature T = 150 K the island shapes are triangular with their sides being 
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Figure 7. Temperature dependence of the edge diffusion rates for atom diffusion along the {100} 
step by hopping with Fp — 2.5 x 10^^ s~^ (sohd Hne), and along the {111} step by exchange with Fq 
= 2.5 X lO^-^ s-i (dash-dotted line). 
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{100} steps. Increasing the temperature to T = 200 K a transition from triangular 
to hexagonal shape occurs and for T = 250 K the islands become triangular again. 
However, at this temperature they are mainly bounded by {111} steps. 

To understand the island shapes in the temperature regime between 150 and 250 K, 
we consider the mobility of the adatoms along the steps (at such temperatures the 
adatoms at the step edges cannot leave the steps): The lower the migration proba- 
bility along a given step edge, the higher the step roughness and faster the speed of 
advancement of this step edge. As a consequence, this step edge shortens as a result 
of the growth kinetics and eventually it may even disappear. Since diffusion along the 
densely packed steps on the (111) surface (the {100} and {111} facets) is faster than 
along steps with any other orientation this criterion explains the presence of islands 
which are mainly bounded by {100} or {111} steps. The same argument can be ex- 
tended to the diffusion along the two close-packed steps and applied to the triangular 
islands at T = 150 K, where the energy barrier for the diffusion along the {111} facet 
is larger and thus the {100} steps survive so that triangular islands with {100} sides 
are obtained. By considering the energy barriers we would expect only these islands, 
until the temperature regime for the thermal equilibrium is reached. However, as noted 
in the introduction, the diffusion of adatoms is not only governed by the energy bar- 
rier but also by the effective attempt frequency. For A1/A1(111) the effective attempt 
frequencies have not been calculated, but the analysis of Ref. 35 proposes that the ex- 
change process should have a larger attempt frequency than the hopping process. The 
results displayed in Fig. ^ are obtained with 1.0 x 10^^ s~^ for the diffusion on the fiat 
surface, 2.5 x 10^^ s~^ for the jump along the {100} step, and 2.5 x 10^^ s~^ for the 
exchange along the {111} step. These effective attempt frequencies are the only input 
of the KMC not calculated explicitly by DFT, but were estimated from the theoretical 
PES as well as from experimental data for other systems. In Fig. |^ the edge diffusion 
rates along the two steps are plotted as a function of the reciprocal temperature. At 
lower temperatures the energy barrier dominates the diffusion rate, but at T = 250 K 
the attempt frequencies start to play a role and lead to faster diffusion along the {111} 
facet than along the {100} one. Thus, the latter steps disappear and only triangles 
with {111} sides are present. The roughly hexagonally shaped islands at T = 200 K 
are a consequence of the equal advancement speed for the two steps at that temper- 
ature. Obviously, the temperature-dependence of the growth shapes found in Fig. ^ 
is crucially determined by the ratio of the two diffusivities and in particular by the 
temperature at which the two lines of Fig. |^ cross. If the difference were only one 
order of magnitude, the crossing would be at a temperature that is too high (namely 
at T = 505 K). The formations of fractals (Fig. ^, upper left) and of {100} step tri- 
angles would still occur. Obviously, the importance of the attempt frequencies should 
receive a better assessment through accurate calculations, and work in this direction 
is in progress. While no experimental data for A1/A1(111) are presently available we 
note that a similar sequences of islands as obtained above has been observed for Pt on 
Pt (111) by Michely et al.}^ 
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